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In many applications, the two-dimensional trajectories of fluid particles are avail- 
able, but little is known about the underlying flow. Oceanic floats are a clear ex- 
ample. To extract quantitative information from such data, one can measure single- 
particle dispersion coefficients, but this only uses one trajectory at a time, so much 
of the information on relative motion is lost. In some circumstances the trajectories 
happen to remain close long enough to measure finite-time Lyapunov exponents, but 
this is rare. We propose to use tools from braid theory and the topology of surface 
mappings to approximate the topological entropy of the underlying flow. The pro- 
cedure uses all the trajectory data and is inherently global. The topological entropy 
is a measure of the entanglement of the trajectories, and converges to zero if they 
are not entangled in a complex manner (for instance, if the trajectories are all in a 
large vortex) . We illustrate the techniques on some simple dynamical systems and on 
float data from the Labrador sea. The method could eventually be used to identify 
Lagrangian coherent structures present in the flow. 

Keywords: topological chaos, dynamical systems, Lagrangian coherent structures 

Consider particles floating on top of a fluid. We can foflow their trajectories, 
either with a camera or by computer simulation. If we then plot their position in 
a three-dimensional graph, with time the vertical coordinate, we get a 'spaghetti 
plot,' which contains information about how entangled the trajectories are. We 
discuss how to measure the level of entanglements in terms of topological en- 
tropy, and the interpretation of the results. This provides a straightforward 
method of estimating the level of chaos present in a system. This approach 
could also be used to determine if some trajectories remain together for a long 
time, and are thus part of a Lagrangian coherent structure. 



I. INTRODUCTION 
A. Floats in the ocean: an example 



Figure 1(a) shows the trajectories of ten floats released in the Labrador sea, for a period 
of a few months. The principal reason to release such floats is the data they measure 
and transmit back - temperature, salinity, pressure, etc. But the actual trajectories of 
the floats are also important, since they tell us something about large-scale transport in 
the ocean, a crucial component in understanding global circulation. From a single float 
one can deduce the single-particle dispersion coefficient, a crude measure of how quickly a 
float wanders away from its release point. However, it is better to measure quantities that 
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FIG. 1: (a) 10 floats from Davis' Labrador sea data [T]. (b) The growth of L(u) for the 10 floats, 
with a fit of the topological entropy. For longer times, the floats leave the Labrador sea and L{u) 



becomes constant. The details of the analysis are presented in Section IV 



involve several floats [2j. For instance, if floats happen to start near each other, then we can 
see how quickly they separate and measure flnite-time Lyapunov exponents [3], which are 
linked to chaotic advection But if the floats are nowhere near each other, then a more 
global quantity is needed. In this paper we propose to examine the 'braid' deflned by the 
trajectories and to measure their degree of entanglement. (All these terms will be deflned 
more precisely.) The number we get out of this is called the braid's topological entropy. 



Figure 1(b) shows a measure of the entanglement of the ten floats as a function of time, 



with an exponential flt: the growth rate is the topological entropy. (For longer times, the 
curve levels off because floats start leaving the Labrador sea, which is not really a closed 
system.) Much like a Lyapunov exponent, the topological entropy gives us a characteristic 
time for the entanglement of the floats in the Labrador sea, here about 1/0.02 ~ 50 days. 



B. Topology and trajectories 



Since the original paper of Boyland et al. [5], topological techniques in fluid dynamics 
have been applied to free point vortices [6J, flxed blinking vortex |71|8], rod stirring devices |9l 
[TOl 111] [T2| IT3] . and spatially-periodic systems [HI |15]. (See [16] for a brief review.) More 
recently, the emphasis has shifted to locating periodic orbits that play an important role in 
stirring [H] — so-called ghost rods — and even to the manufacture of such orbits [T7] . 
Most of these papers study periodic motions of rods or particle orbits. For many practical 
applications, however, periodic motion is not directly observable, since most such orbits 
are highly unstable. Hence, some authors have examined random braids [3, Ull HE] [TH] [T9| 
composed of arbitrary chaotic trajectories. (There is also related literature from the knot 
theory perspective — see for example [20].) 

The goal of the present paper is to give concrete techniques that can be used to obtain 
topological information from particle trajectories. The mathematical details are glossed over: 
the emphasis is on usability. Implementation details are discussed, and some sample Matlab 
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programs are presented in an appendix. The hope is that this will make these techniques 
more accessible to those with little or no background in braid theory and topology. 

The principal measurement we extract from a braid is its topological entropy. This en- 
tropy is closely related to the traditional Lyapunov exponent, except that being a topological 
quantity it is not sensitive to the size of the sets on which chaos is occurring. This is both 
a weakness and a strength: it does not tell us everything we might like to know, but on 
the other hand the topological entropy is easy to compute from crude data. This is in con- 
trast to Lyapunov exponents, which require at the very least detailed knowledge of particle 
trajectories that start close together, and at best the velocity field and its gradient. When 
dealing with, for example, data from oceanic floats (as we will later in this paper), being 
able to compute a Lyapunov exponent is a rarity. 

There is also a philosophical point that bears some discussion. The viewpoint of the 
present paper is that given particle trajectory data, a useful thing to quantify is how 'entan- 
gled' the particle trajectories are. This can be done from the particle data directly, without 
worrying about the underlying flow. By contrast, Lyapunov exponents are defined locally 
and are sensitive to the smooth structure of the flow. It is exactly the (presumed) smooth 
nature of the flow that connects local information to a global quantity such as the Lyapunov 
exponent, but one is left wondering why we should care about the local picture at all in 
practical situations. The topological viewpoint presented here is an attempt to sidestep this 
and focus directly on global information. 

We begin in Section [IT] by a short introduction to braid theory, surface dynamics, and 



their connection to dynamical systems. In Section II B we show how to extract braids from 
particle trajectory data. Section III is devoted to topological entropy: in Section III A we 



discuss its connection to flows, and in Section III B we show how to measure it from a braid 



(for a braid corresponding to periodic orbits). In Section IV we introduce random braids, 
and again show how to measure entropies. As an application, we calculate the entropy for 
floats in the Labrador sea, as presented in Fig. [T] We offer some concluding remarks in 
Section |Vl 



II. BRAIDS 
A. Physical and Algebraic braids 



First we describe intuitively how braids arise. Figure 2(a) shows the orbits of 4 particles 



in a circular two-dimensional domain. The particles might be fluid elements, solutions of 



ordinary differential equations, or physical particles at the surface of a fluid. Figure 2(b) 



shows the 'world line' of the same orbits: they are plotted in a three-dimensional graph. 



with time flowing vertically upwards. The diagram in Fig. 2(b) depicts a physical braid, 
made up of four strands. No strand can go through another strand as a consequence of the 
deterministic motion of the particles (they never occupy the same point at the same time). 
Moreover, the mathematical definition of a braid requires that strands cannot 'loop back': 
here this simply means that the particles cannot travel back in time. We will say that two 
braids are equivalent if they can be deformed into each other with no strand crossing other 
strands or boundaries. Throughout this paper, we will be interested in characterizing the 
level of 'entanglement' of trajectories. 

Since we can move the strands, is convenient to draw braids in a normalized form, as 



shown in Fig. 2(c) for the braid in Fig. 2(b) Such a picture is called a braid diagram. The 



(a) 
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FIG. 2: (a) The orbits of four particles in a circular two-dimensional domain; (b) The same orbits, 
lifted to a space-time diagram in three dimensions, with time flowing from bottom to top. (c) The 
standard braid diagram corresponding to (b). 



important thing is that we record when crossings occur, and which particle was behind and 
which was in front. It matters little how we define 'behind,' as long as we are consistent 



(see Section II B for practical considerations). In Fig. 2(c) the horizontal dashed lines also 
suggest that we can divide the braid into a sequence of elementary crossings, known as 
generators. Figure 3(a) shows the definition of cTj, which denotes the clockwise interchange 
of the ith and {i + l)th strands, keeping all other strands fixed. Note that the index i is the 
position of the strand from left to right, not a label for the particular strand. For n strands, 
we have n — 1 distinct generators. 

Figure 3(a) also shows the counterclockwise interchange of two strands, denoted by the 
operation a^^ . The justification of the 'inverse' notation is evident in Fig. 3(b) if we 
concatenate a,: and a"^, then after pulling tight on the strands we find that they are dis- 



entangled. We call the braid on the right in Fig. 3(b)| the identity braid. In fact, the set 
of all braids on a given number n of strands forms a group in the mathematical sense: the 
group operation is given by concatenation of strands, the inverse by reversing the order and 
direction of crossings, the identity is as described above, and it is clear that concatenation is 
associative. This group is called the braid group on n strands, also known as the Artin 
braid group. 

The braid group Bn is generated by the set {di, . . . ,an-i}. this means that any braid 
in Bn can be written as a product (concatenation) of cii's and their inverses. The braid group 
is finitely- generated, even though it is itself infinite: only a finite number of generators give 
the whole group. To see that the braid group contains an infinite number of braids, simply 
consider af, for k an arbitrary integer: no matter how large k gets, we always get a new 
braid out of this, consisting of increasingly twisted first and second strands. 



We have now passed from physical br aids, as depicted in Fig. 2(b) , to algebraic braids. 
The algebraic braid corresponding to Fig. 2(c) is a^^(T2^(r^^(r2<7i, where we read generators 
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FIG. 3: (a) Braid group generator ai, corresponding to the clockwise interchange of the string at 
the ith and {i + l)th position, counted from left to right. Its inverse a^^ involves their counter- 
clockwise exchange, (b) The concatenation of Cj and a^^ gives the identity braid. 



from left to right in time (beware: conventions differ). In essence, an algebraic braid is simply 
a sequence of generators, which may or may not come from a physical braid. How can we 
guarantee that physical braids and algebraic braids describe the same group? We need to 
be mindful of relations amongst the generators that arise because of physical constraints. 
For example. Fig. 4(a) shows a relation amongst adjacent triplets of strands. Staring at 



the picture long enough, and allowing for the deformation of strands without crossing, the 
reader can perhaps see the that braids in Fig. 4(a) are indeed equal. Hence, the algebraic 
sequence (jjCTj+iaj must be equal to (Ji+io-jO-j+i, if the generators are to correspond to physical 



braids. Another, more intuitive relation is shown in Fig. 4(b) generators commute if they 
do not share a strand. In summary, we have the relations 



if \i- j\ = 1, 
if > 1, 



(la) 
(lb) 



amongst the generators. Artin [21] proved the surprising fact that there are no other relations 
satisfied by the generators ai, except for those than can be derived from ([T| by basic group 
operations (multiplication, inversion, etc.). The generators {ai, . . . , cr„_i} together with the 
relations ([T| define the algebraic braid group, which we also denote With these relations, 
the groups of physical and algebraic braids are isomorphic. 

A consequence of the relations ([T]) is that it may not be immediately obvious that two 
algebraic braids are equal. For instance, the braids aia2 and a2a2(Tia2(7i^ cti^ are equal, 
since 



(Ti(T2 



i(Tl(T2(Ti (T 



cr2(cricr2(Tij(Ti 



a 2(J 2(^1(^2(^1 ^^1 ^ 



This 'braid equality' problem has seen many refinements: the original solution of Artin [21] 
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FIG. 4: Braid group relations (see Eq. ([T|): (a) relation for three adjacent strands; (b) commutation 
relation for generators that don't share a strand. 



has computational complexity exponential in the number of generators, but modern tech- 
niques can determine equality in a time quadratic in the braid length |22l ESI El] . 

B. Extracting the braid from a flow 

The first step in obtaining useful topological information from particle trajectories is to 
compute their associated braid, essentially going from the physical picture in Fig. 2(b)| to 



the algebraic picture in Fig. |2(c) A simple method to do this was originally described in [7] , 
but is also implicit in earlier work such as [HI [25] (see also [26] for a related technique). 

We start with trajectory information for n particles over some time. We first project the 
position of the particles onto any fixed projection line (which we choose to be the horizontal 
axis), and label the particles by i = 1,2, ... ,n in increasing order of their projection. A 
crossing occurs whenever two particles interchange position on the projection line. A crossing 
can occur as an "over" or "under" braid, which for us means a clockwise or counterclockwise 
interchange. These interchanges correspond to the braid group generators introduced in 
Section EAl 

Assuming a crossing has occurred between the ith and {i + l)th particles, we need to 
determine if the corresponding braid generator is cTj or 0"^"^. We look at the projection of 
the iih and {i + l)th particles in the direction perpendicular to the projection line (the 
vertical axis in our case). If the iih particle is above the (i + l)th at the time of crossing, 
then the interchange involves the group generator (Xj (we define "above" as having a greater 
value of projection along the perpendicular direction). Conversely, if the ith. particle is below 
the {i + l)th at the time of crossing, then the interchange involves the group generator . 
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FIG. 5: Detecting crossings: (a) Two possible particle paths that are associated with different 
braid group generators; (b) Two crossings that yield no net braiding. The projection line used to 
detect crossings is shown dotted at the bottom, and the perpendicular lines used to determine the 
braid generator are shown dashed. 



Figure 5(a) depicts these two situations. 



The net result is that from the n particle trajectories we obtain a time-ordered sequence 



of the generators ai and 



n — 1. We call this sequence the braid of the 



trajectories. We also record the times at which crossings occur, so each generator in the 
sequence has a time associated with it. 

Remarks: 

• The method just described might seem to detect spurious crossings if two well- 
separated particles just happen to interchange position on the projection line several 
consecutive times in a row, as shown in Figure 5(b) However, this would imply a 
sequence of ai and a^^ braid generators, since which particle is the ith one changes at 
each crossing. When composed together the generators for these crossings cancel. 



We give a simple Matlab implementation of the method in Appendix A 1 The pro- 
gram gencross detects crossings in trajectory data; it makes an effort to resolve 
multiple simultaneous crossings (up to triple crossings), but will complain if it gets 
confused. More sophisticated code can be written that re-interpolates the trajectory 
as needed to detect crossings. 

When the system has symmetries, such as when several periodic orbits lie on the same 
line, there are 'bad' choices of projection line where it is impossible to resolve the order 
of crossings, since orbits cross at exactly the same time. Displacing the projection line 
a little cures this. 

If the braid is truly periodic, that is, if all the particles return to their initial configu- 
ration, then changing the projection line changes the braid, but only by conjugation, 



which does not affect the entropy [6] (Section III). 
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If the trajectories are not periodic, then the method does not define a braid in the 
traditional sense where all the strands return to the same initial configuration. This 
is inconsequential to our purposes: all that matters is the order along the projection 
line (see also [2S])- The choice of projection line changes the braid beyond simple 
conjugation, but this only creates an error in a small, finite number of generators, 
which is not important when considering long braids and does not asymptotically 



affect the entropy (Section III) 



If the braid is generated from chaotic trajectories, then missing a few crossings (due 
to, say, gaps in the data) is fine as long as the trajectories are long enough. 



III. TOPOLOGICAL ENTROPY 

In Section |TT] we described how a set of trajectories in the two-dimensional dynamical 
system can be described as a braid in a three-dimensional space-time diagram. In this 
section we will describe further how this braid relates to topological information for the 
underlying fiow. 

It is worth noting that braids are not always interpreted in terms of trajectories: they 
arose first and are still studied as independent geometrical and algebraic objects. The rea- 
son they take center stage in the present study is through their connection to mappings of 
surfaces (mapping class groups) . The Thurston-Nielsen theory [271 ESI [2HI EDI EI] classifies 
mapping of surfaces according to whether they can be "deformed" to each other in a topo- 
logical sense. Braids provide a convenient way of labeling the isotopy classes that result. So 
even though we will often speak here of the braid as being the primary object of interest, 
we are really using techniques that apply to the class of mappings labeled by a braid. 



A. Entropy of a flow 

Ultimately, we want to measure the topological entropy of a system directly from a braid 
of trajectories. Before we do this, we discuss the meaning of topological entropy of a fiow or 
map. The topological entropy of a dynamical system measures the loss of information under 
the dynamics. It is closely related to the Lyapunov exponent, which measures the time- 
asymptotic rate of separation of neighboring trajectories. But it is in some sense a cruder 
quantity, since it does not require a notion of distance. A positive entropy is associated 
with chaos, though it tells us nothing about the size of the chaotic region. The topological 
entropy is an upper bound on the largest Lyapunov exponent of a fiow. The two are equal 
only for very simple systems where stretching is uniform. 

Though there are more fundamental ways to define it, we shall take our working definition 
of entropy to be the asymptotic growth rate of material lines [32]. It is fairly straightforward 
to measure this numerically, given a sufficiently accurate velocity field. We simply choose 
an initial material line and follow it for some time, interpolating new points as the line gets 
longer. Figure |6] shows such a line for a numerical simulation of a stirred viscous fiow. Note 
how the exponential growth rate is very sharply defined. The topological entropy is the 
supremum of the growth rate over all such loops, but in practice almost any nontrivial loop 
{i.e., that spans the domain) will grow exponentially at a rate haow 

In practical applications we often do not have access to an accurate representation of 
the velocity field. This is where braids come in, as a way of approximating the topological 
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FIG. 6: The length of a material line as it is advected by a flow for many periods. The exponential 
growth rate is very well defined (fitted line), even though there are identical small oscillations at 
each period. 



entropy. As we will see in Section IIIB the braid provides a lower bound on the flow's 
topological entropy. 



B. Entropy of a braid 

Figure [7] illustrates how the motion of n = 3 point particles can be used to put a lower 
bound on the topological entropy, defined here as the growth rate of material lines or loops. 
Here, the point particles undergo a motion described by the braid aiO"^^. Two full periods 
are shown. Notice that an initial loop that is 'caught' on the particles is forced to follow 
along, since determinism implies that it cannot occupy the same point in phase space as 
the particles. In fact, a straightforward calculation [5] shows that for this braid the total 
length of the loop must grow exponentially at least at a rate log((3 + -\/5)/2) per period. 
We call this rate the topological entropy of the braid, h^l^^^, to distinguish it from the true 
topological entropy of the flow, h^ow, as defined in Section [III A[ We have 



^flow > ^'■"^ 



braid 



(2) 



for any braid obtained from the motion of n particles in the flow. Typically, the more 



particles are included in the braid, the closer h^l^^^ is to /iflow [LM- Note that /i^riid ^braid 
are always zero. 

An essential property of h^}^^^ is that the growth rate of the loop is independent of specific 
details: for instance, if the particles are not equally spaced, or if the loop is 'tightened' around 
the particles, then the length will change, but the asymptotic growth rate will not, because 
all these changes amount to an additive constant in the logarithm, which gets divided by a 
large time. 



,(1) 



,(2) 
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FIG. 7: For n = 3 particles, a loop is initially wrapped around the second and third parti- 
cles. The generator cJi is applied, interchanging clockwise the first and second particles, followed 
by cr^^, which interchanges the second and third particles counterclockwise. The loop is forced to 
stretch under this application. The final loop on the lower right underwent two applications of the 
braid uicr^^. If we keep applying this braid, the length of the loop will grow exponentially. 




FIG. 8: A non-intersecting closed loop wrapped around n = 5 particles. Up to trivial deformations, 
the loop can be reconstructed by counting intersections with the vertical lines. In terms of the 
crossing numbers defined in Fig. (9) this loop has z^i = z^2 = 2, 1^3 = 1^4 = 4, /ii = 2, ^3 = 1, ^us = 4, 
A*2 = = 0, 114 = 3, and Dynnikov coordinate vector u = (—1, 1, —2,0, —1,0) (see Eq. Q). 

We are now faced with a task: given a sequence of generators (Xj, measured in some way or 
obtained numerically from a flow, what is h^l^^^? The method used in [7], based on a matrix 
representation of the braid group, only provides a lower bound on the braid entropy. An 
accurate and efficient computation has since become a lot simpler due to a new algorithm 
by Moussafir |33], who uses a set of coordinates to encode a loop. We describe this briefly 
below; for more details see [331 EH [35]. The reader who is mostly interested in using the 
method can skip to the end of the section to Procedure [l| 

The basic idea is simple: consider the closed loop in Fig.[8| which is wrapped around n = 5 
particles. The loop does not intersect itself, so in two dimensions the allowable paths it 
can follow around the particles are far from arbitrary. The amazing fact is that we can 
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FIG. 9: Definition of tlie crossing numbers fii and u^. The for i odd count crossings above a 
particle, and below a particle for i even. The z/j count crossings between particles. 



reconstruct the entire loop, or at least the way it is threaded around the particles, by 
counting how many times it intersects the vertical lines in Fig. |8j 

In Fig. |9] we give specific labels to the crossing numbers. For n particles, /i2i-3 (odd 
index) gives the number of crossings of a loop above the ith particle, and ^2i-2 (even index) 
below the same particle. The number z/j counts the crossings between particle i and i + 1. 
We have a total of 3n — 5 crossing numbers. 

This set of crossing numbers (which are all non-negative) can be reduced further: define 

0-1 = \ (yU2i - yU2i-i) , h = \{yi- z/j+i) (3) 
for z = 1, . . . , n — 2. The vector of length {2n — 4), 

u = (ai, . . . ,a„_2,6i, . . . ,6„-2) (4) 

is called the Dynnikov coordinates of a loop. (As far as we know, this specific encoding was 
originally introduced by Dynnikov [31] , but it is implicit in earlier work of Thurston, Dehn, 
and others.) The components and hi are signed integers. They can be used to exactly 
reconstruct the loop but we shall not need to do this here. This set of coordinates is 
minimal: it is not possible to achieve the same reconstruction with fewer numbers. 

We can also obtain the minimum number of intersections L{u) of the loop with the 
horizontal line through the particles |33j : 

n— 3 n— 1 

L{u) = \ai\ + |a„_2| + - ai\ + ^\bi\ , (5) 

1=1 i=Q 

where bo and 6„_i can be obtained from the other coordinates as [35j 

n-2 

bn-i = -bo-^bi. (6) 
i=i 



i-i 

6o = — , max ( \ai\ + bf + bj 



l<i<n~2 



12 



The formula for L{u) is encoded in the Matlab function loopinter in Appendix A 2 For 



example, the loop in Fig. |8] intersects the horizontal axis (the line through all the particles) 
12 times. The crucial observation, which will allow a simple computation of h^}^^^, is that if 
the length of the loop grow exponentially, then L{u) also grows exponentially at the same 
rate [271 E3]. 

Now that we have a way of encoding any loop, we need to find how the loop is transformed 
by a braid. What makes all this work is that there is a very efficient way of doing this: given 
a loop encoded by u as in Eq. Q, each generator of the braid group ai simply transforms 
these coordinates in a predetermined manner. (Mathematically, this defines an action of the 
braid group on the set of Dynnikov coordinates.) We call these transformations the update 
rules for a generator. 



The update rules are straightforward to code on a computer. (See Appendix A 2 for 
a Matlab implementation.) To express them succinctly,^ first define for a quantity / the 
operators 

/+:=max(/,0), /" := min(/, 0). (7) 

After we define 

Q_i = ai_i -tti- bf + , (8) 
we can express the update rules for the cTj acting on u = (ai, . . . , a„_2, 6i, • • • , bn-2) as 



a' 



- ai_i - 6+1 - {bf + , (9a) 

K-i = h + cr , , (9b) 

a'i = ai- br - (for^ - Ci_i)~ , (9c) 

= k-i - c- 1 , (9d) 

for i = 2, . . . , n—2. For this and the following update rules, all the other unlisted components 

of u are unchanged under the action of Uj or a~^. The leftmost {i = 1) and rightmost 
{i = n ~ 1) generators require special treatment, having update rules 

a'l = -&! + (ai + 6+)+, (10a) 

b[ = ai + bt . (10b) 



for (Ji, and 



I 



-bn-2 + (an-2 + b^_^ , (11a) 
b'n_2 = an~2 + &;_2 • (lib) 



for 

We need to give separate update rules for the generators a~^. With the definition 

di_i = ai_i - + bf - br_^ , (12) 



^ We are using the numbering scheme of [3S], but the notation of Also, we define generators ai as 
clockwise interchanges rather than counterclockwise. 
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the update rules for the ^ are 



= ctj-i + + - rfj-i) , (13a) 

=h- dU , (13b) 

a'. = ai + + [hj_^ + di-i) , (13c) 

b[ = + dU , (13d) 



for z = 2, . . . , 77, — 2. We also have 



% = bi — (b^ — ai)^ , (14a) 
b[ = 6+ - ai . (14b) 



for cr^ \ and 



a. 



n-2 — " ^ a-n-2) , (I5a) 

= K-2 - an-2 • (I5b) 



for a^\. 

Update rules of this form are known as piecewise-linear: once the mins and maxes are 
resolved, what is left is simply a linear operation. However, the mins and maxes are what 
keeps this from being a simple linear algebra problem and make the braid dynamics so rich. 

Here then is a recipe for computing /i[,"aid; ^be topological entropy of a braid of n particle 
trajectories: 



Procedure 1 (Entropy of periodic braid) 

1. Start with an arbitrary initial loop, encoded as a vector u (Eq. Set m to 0; 



2. For each generator in the braid, use the appropriate update rule ([8j)-(15) to modify u; 

3. Compute the intersection number L{u) using Eq. ([s]); 
4- Repeat steps^^^for all generators in the braid; 



5. Add 1 to m; Calculate hllll-, = m ^ log L(u); 



6. Repeat steps to until converges in step 



Remarks: 



The procedure above assumes that the braid is periodic, i.e., is obtained from periodic 
orbits of the flow. In Section HV] we will discuss how the method differs for random 
braids obtained from sampling arbitrary trajectories (which don't necessarily repeat). 

The dimension of /i[,"^id inverse time, where the unit of time is the period over which 
the braid is repeated. 

Even though the discussion so far has described -u as a vector of integers, the initial 
condition for u in step [T] can in practice be chosen to be a random set of real numbers. 
(This called a projectivized version of the coordinates [35].) 
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FIG. 10: For the braid a^^a2^a^^a2(Ti in Fig. [2] (a) Entropy estimate log L{u) as a function 
of period m using Procedure [!} (b) Error (deviation from the true entropy ~ 0.83144), showing 
the 1/m convergence (dashed hne). 



As is typical of such exponential growth calculation, it is possible that the compo- 
nents of u becomes so large that they overflow double-precision arithmetic. In that 
case standards 'renormalization' techniques can be used: divide -u by a large con- 
stant I/overflow5 but keep track of how many times this division was done. Then add 
that multiple of log -^overflow to the logarithm in step [5] Another option is to use real 
or integer arbitrary precision arithmetic, but this slows down the calculation. 

We stress that Moussafir's technique for the computation of a braid's entropy is ex- 
tremely rapid compared to previous methods, which typically use train tracks and the 
Bestvina-Handel algorithm [36J, or combinatorial methods |2S]- The rapidity arises 
from the fact that the algorithm keeps a bare minimum of information (the vector u) 
to express the topology of an arbitrarily long curve. The Bestvina-Handel algorithm, 
however, gives more information about the braid (such as the existence of invariant 
curves - see Section Ivl) . 



The speed of convergence of this procedure is discussed in [HI [33] . As an example. 
Fig. |10(a) shows the result of applying the procedure to the braid in Fig. [2| and Fig. 10(b) 



shows the convergence rate to the exact entropy. 



IV. RANDOM BRAIDS 

From the point of view of data analysis, looking at periodic braid is not general enough. 
Most periodic orbits in a dynamical system are unstable, and thus they cannot be detected 
directly. The trajectories we have access to are typically chaotic. Nevertheless, we can ask 
what the braid corresponding to a set of orbits tells us about a dynamical system. The 
answer is that its entropy approximates the 'true' topological entropy of the flow, and the 
approximation gets better as more particles are added. 
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There are two ways to analyze random braids generated by chaotic trajectories: without 
and with ensemble averaging. 'Without averaging' means that we have a single realization to 
study, say n trajectories integrated or measured up to some final time. Unless the final time 
is extremely long, this is not very accurate. 'With averaging' means that we have the luxury 
of repeating the experiment several times, following the same number of trajectories at each 
realization (assuming the fiow is the same for each realization, at least in a statistical sense). 
We then average over the total number of realizations of the experiment, in the manner 
described below. 



A. Entropy without averaging 

Let us first describe the procedure without averaging: we assume that we have obtained 
a sequence of generators from the trajectories of n particles, as well as the time at which 



each crossing occurs (see Section II B). In the examples presented here, those trajectories 
were either computed from randomly-selected initial conditions, or they were obtained from 
measured data (the oceanic fioats). 



Procedure 2 (Entropy of random braid, without averaging) 

1. Start with an arbitrary initial loop, encoded as a vector u (Eq. 



2. For each generator in the braid, use the appropriate update rule ([8j)-(15) to modify u; 

3. Compute the intersection number L{u) using Eq. ([s]); 

4- Plot log L{u) versus t, where t is a vector of times when each crossing occurs; 

5. Repeat steps to until we can fit a line in step 



Remarks: 

• Since the braid is random, we must keep track of the time when crossings occur. 

• Most of the remarks from the periodic Procedure [l] of Section III B still apply. 



Figure [Tl] shows an example of applying Proce dure [2 to n = 3 particles advected by a blinking 



vortex fiow [U [71 [H] in the regular regime (Fig 11(a) ) and chaotic regime (Fig. 11(b) ). In the 



first case, the growth of L{u) is roughly linear, so the entropy is zero. In the second case the 



growth is exponential. Note that the integration time is quite long, and in Fig. 11(b) L{u 
becomes enormous. For such long integration time, the fit for the entropy is good. 



B. Entropy with averaging 

To get a more accurate measurement of /i[,"aid random braids, ensemble averaging is 
desirable, if we have that luxury. To implement this, we integrate a set of n trajectories A^reai 
times, randomizing the initial condition for each realization. We obtain a list of n x A^reai 
trajectories for times < t < tmax, from which we compute A^reai braids and vectors of 
crossing times. To do the averaging, we need to be able to compare L{u) at the same 
times for each braid, but since crossings occur at different times we cannot do this directly. 
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FIG. 11: (a) The number of crossings of a loop, versus time, for the blinking vortex flow in the 
regular regime (circulation T = 0.5, co-rotating vortices; see [7j). The growth of L(u) is linear, (b) 
Same as in (a), but in the chaotic regime (circulation T = 16.5, counter-rotating vortices; see [7]). 
The vertical axis is now on a log scale; the slope of the line gives the braid entropy (Procedure |2]). 



We instead break up the total time interval into equal subintervals of length At, and for 
each subinterval and each realization we record L{u) up to time qAt, where q is an integer 
with < qAt < tmax- We finally obtain A^reai lists of [t^ax/At] (square brackets denote the 
integer part) intersection numbers L{u), all sampled at the same times corresponding to 
each subinterval. The whole procedure is summarized as: 



Procedure 3 (Entropy of random braid, with averaging) 

1. Start with iVreai lists of intersection numbers L{u) and their crossing times, generated 
following Procedure [1[ steps 

2. For each realization, record the intersection numbers up to fixed times qAt, 

< qAt < tmax/ 

3. At each time qAt, compute the average (logL(it)) over all A^reai realizations; 
4- Plot {log L{u)) versus qAt, < qAt < tmax; and fit a line to get /i[,"aid- 



Remarks: 



We average logL('u) rather than L{u): not only does this ensure that Procedures [2] 
and [3] give the same entropy, but it also leads to smaller fluctuations. 

For best results, the number of subintervals [tmax/ At] has to be large enough to get a 
good fit, but small enough that there are several crossings within each subinterval of 
length At. 



Figure IV B shows an example of applying Procedure |3] with A^reai = 50 realization of n = 3 
particles advected by the same blinking vortex flow as for Fig. 11(b), Notice that the fit 
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FIG. 12: Similar plot to Fig. 11(b), but after averaging (logL(it)) over A'reai = 50 shorter trajec- 
tories (tmax = 100 time units rather than 1000). The average is plotted at each At = 10 time units 
(see Procedure [s]) . The individual trajectories are shown in the background. 

is much better, even though the integration time is shorter. We used [tmax/c^^] = 10 time 
subintervals of length At = 10. An expUcit example in Matlab (for the Duffing oscillator) 



is given in Appendix A 3 



Oceanic floats 



As a more practical application of random braids, we consider data for oceanic floats 



in the Labrador sea (North Atlantic) [I], discussed in the introduction (Section I A). The 
position of ten floats for a few months is shown in Fig. 1(a) Note that the float trajectories 



seem more entangled while they are confined to the Labrador sea (between Greenland and 
Labrador), and some eventually escape. To compute the braid, we linearly interpolate the 
float positions to determine when crossings occur between the n = 10 floats. We then use 
Procedure |2] to compute the entropy, as shown in Fig. 1(b), since ensemble averaging is 



not available here (we only have data for one experiment). We see in Fig. 1(b) that L(u) 



has a convincing exponential regime for about 150 days, after which floats tend to escape 
the Labrador sea and L{u) reaches a plateau. The entropy gives us a timescale for the 
entanglement of floats in the Labrador sea, here about 1/0.02 ~ 50 days. This number is 
easy to obtain from the raw data: there is no need for a model of the velocity field. However, 
the trajectories need to be long enough for a significant number of crossings to occur, and 
localized enough for particles to actually braid. 

More context will be needed to fully understand what it means to say that the timescale 
for entanglement is 50 days. For instance, the method could be benchmarked by following 
tracers in simulations of flows comparable to the Labrador sea, in which braids can be easily 
computed. The measure is also useful for comparing different regions of the ocean. 

Note that there has been previous work on bounding the topological entropy of experi- 
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mental data. Amon and Lefranc [37] have obtained lower bounds on entropy and evidence 
of chaotic behavior in a nonstationary optical system. See also the review by Gilmore |38j . 



V. DISCUSSION 

There are two ways to interpret the data obtained from braids of particles. The first 
is to accumulate data for enough particle trajectories that a good approximation to the 
topological entropy /ifiow is obtained. The drawback to this is that the convergence of h^^^^ 
to /iflow as n gets larger appears to be fairly slow [T^, though more work is needed to 
determine this convergence rate. In this interpretation, the braid approach is seen as a 
practical way of measuring /iaow This interpretation is in the same spirit as a Lyapunov 
exponent. Its main advantage is that a single number is easy to comprehend and compare; 
its main drawback is that a single number doesn't capture the subtleties of a particular 
system. 

The second interpretation is to regard h^l^^^ as the 'n-particle braiding time,' in a similar 

fashion that n-particle correlation functions are measured. Thus, the behavior of h^]^^^ 
with n carries real information, as it tells us the typical time for 3 particle trajectories to 
become entangled, then 4 particles, etc. We might call this the 'spectrum of braid entropies' 
for a dynamical system. The drawback of this approach is that it requires a more careful 
analysis of the data. 

The method presented in this paper is limited to two-dimensional flows. Indeed, a four- 
dimensional braid of three-dimensional particles trajectories is not very useful, as this and 
all higher-dimensional braid groups are trivial (strands can always be disentangled without 
crossing [31]). The best alternative is to lift the trajectories of material lines to sheets in 
four dimensions, but this presents some daunting visualization challenges, and there is little 
developed theory (but see [39]). 

Finally, the topological entropy h^]^^^ is only the crudest piece of information that can be 
extracted from a braid. Many other types of invariants can be computed [IQ]. However, those 
invariants don't necessarily have a clear interpretation in terms of dynamics. A promising 
avenue for obtaining much more precise information on a dynamical system is to find the 



isotopy class of the random braid (see Section III). This would tell us, for instance, whether 



some floats merely orbit each other and thus behave as one 'trajectory' from the point of 
view of braiding and entropy. This is akin to making a braid out of thick rope: even though 
each rope is made up of tiny strands, these contribute to the braid as one large strand. 
This sort of approach could help to identify Lagrangian coherent structures from particle 
trajectory data, by looking for decomposable braids. However, the tools available to do this, 
such as the Bestvina-Handel algorithm [36], are still slow and difficult to use on large braids. 
A promising approach was recently used in [35], but needs to be developed further. 
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APPENDIX A: MATLAB EXAMPLE PROGRAMS 

These Matlab programs can be obtained by downloading the source of this paper at 



http : //arXiv . org/abs/0906 . 364?! 



1. gencross and interpcross 



The function gencross computes the generators and crossing times for particle trajec- 
tories (see Section II B). The function interpcross is a helper function to gencross that 
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interpolates crossings. Both these functions arc simple implementations with few bells and 
whistles: gencross deals with two or three adjacent particles crossing between successive 
timesteps, but it does not attempt to refine the trajectory (by interpolation or integration) 
to resolve crossings. If it gets confused because too many crossings are occurring between 
two successive timesteps, there is not other option but to refine the data further. 



function [varargout] = gencross (t,X,Y) 

■/.GENCROSS Find braid generators from crossings of trajectories. 

■/. G = GENCROSS (T,X,Y) finds the braid group generators associated with 

crossings of particle trajectories. Here T is a column vector of times, 
and X and Y are coordinates of particles at those times. X and Y have 

'/, the same number of rows as T, and N columns, where N is the number of 
particles. A projection on the X axis is used to define crossings. 

■/. 

■/. [G,TC] = GENCROSS (T,X,Y) also returns a vector of times TC when the 
% crossings occurred. 

'/ Find the permutation at each time. 
[Xperm, Iperm] = sort (X, 2); 

dperm = diff (Iperm, 1) ; '/, Crossings occur when the permutation changes, 

icr = f ind(any(dperm,2)) ; % Index of crossings, 
gen = □ ; tcr = □ ; 

for i = 1 : length(icr) 

Order (from left to right) of particles involved in crossing, 
igen = f ind(dperm(icr (i) , : ) ) ; 

j = 1; 

while j < length(igen) 

if "sumCdpermCicr (i) ,igen(j : j+1))) 
•/. 

'/. Crossing involves a pair of particles. 
•/. 

p = Iperm(icr (i) , igen(j : j+1) ) ; % The two particles involved in crossing. 
[tt,dY] = interpcross(t,X,Y,icr(i) ,p(l) ,p(2)) ; 
tcr = [tcr;tt] ; gen = [gen; igen(j)*dY] ; 

j = j+2; 

elseif "■sum(dperm(icr(i) ,igen(j : j+2))) 
•/. 

■/, Crossing involves a triplet of particles. 

■/, Two cases are possible : 

•/. 

if Iperm(icr(i) ,igen(j)) == Iperm(icr(i)+l,igen(j)+l) 
■/. Case 1: ABC -> CAB 

Particles B&C cross first 
p = Iperm(icr(i) ,igen( [j+1 j+2])); 
[tt,dY] = interpcross(t,X,Y,icr(i) ,p(l) ,p(2)) ; 
tcr = [tcr;tt] ; gen = [gen; igen( j+1) *dY] ; 

'/, Particles A&C cross second 
p = Iperm(icr(i),igen([j j+2])); 
[tt,dY] = interpcross(t,X,Y,icr(i) ,p(l) ,p(2)) ; 
tcr = [tcr;tt]; gen = [gen; igen(j)*dY]; 
elseif Iperm(icr(i) ,igen(j)) == Iperm(icr(i)+l,igen(j)+2) 
•/, Case 2: ABC -> BCA 

7, Particles A&B cross first 

p = Iperm(icr(i),igen([j j + 1])); 

[tt,dY] = interpcross(t,X,Y,icr(i) ,p(l) ,p(2)) ; 

tcr = [tcr;tt] ; gen = [gen; igen(j)*dY] ; 

Pairticles A&C cross second 
p = Iperm ( icr (i), igen ([j j+2])); 
[tt.dY] = interpcross(t,X,Y,icr(i) ,p(l) ,p(2)) ; 
tcr = [tcr;tt] ; gen = [gen; igen( j+1) *dY] ; 
else 



22 



error (' something' ' s wrong with triple crossing — increase resolution') 
end 

j = j+3; 
else 

error ('too many simultaneous crossings -- increase resolution') 
end 
end 
end 

varargout{l} = gen; 

if nargout > 1, varargout{2} = tcr; end 



function [tc.dY] = interpcross (t ,X , Y, itc ,pl ,p2) 
y.lNTERPCRQSS Interpolate a crossing. 

7. [TC.DY] = INTERPCROSS (T,X,Y, ITC, PI, P2) is a helper function for 

y, GENCRQSS. The input is the data T,X,Y (described in the help for 

y, GENCRQSS); the index ITC of the time of crossing (i.e., the particles 

y, cross between T(ITC) and T(ITC+1); and the indices PI and P2 of the two 

y particles that are crossing. INTERPCROSS returns the interpolated 

y crossing time TC, as well as DY (the sign of the difference in Y 

y coordinates) which determines the sign of the generator. 

y Refine crossing time and position (linear interpolation) . 

dt = t(itc+l) - t(itc); Time interval. 

y Particle velocities in that interval. 

Ul = (X(itc+l,pl) - X(itc,pl)) / dt; 

VI = (Y(itc+l,pl) - Y(itc,pl)) / dt; 

U2 = (X(itc+l,p2) - X(itc,p2)) / dt; 

V2 = (Y(itc+l,p2) - Y(itc,p2)) / dt; 

y Interpolated crossing time and Y coordinates. 

dtc = -(X(itc,p2) - X(itc,pl)) / (U2-U1) ; 

tc = t(itc) + dtc; 

Ylc = Y(itc,pl) + dtc*Vl; 

Y2c = Y(itc,p2) + dtc*V2; 

dY = sign(Ylc-Y2c) ; 

y The sign of Ylc-Y2c determines if the crossing is g or g"-l. 
if dY == 0, 

error (' can' 't resolve sign of generator — increase resolution'); 
end 



2. loopsigma and loopinter 



tion 



The function loopsigma applies a sequence of braid group generators to a loop (Sec- 
Eqs. ([8|-([l5))). The function loopinter computes L{u), the minimum number 



IIIB 



of intersections of a loop with the horizontal axis (Section [mB| Eq. (|5 



function up = loopsignia(ii ,u) 

yLOOPSIGMA Act on a loop with a braid group generator sigma. 
y UP = LOOPSIGMA (J,U) acts on the loop U (encoded in Dynnikov coordinates) 
y with the braid generator sigma_J, and returns the new loop UP. J can be 
y a positive or negative integer (inverse generator) , and can be specified 
y as a vector, in which case all the generators are applied to the loop 
y sequentially from left to right. 



n = length (u)/2 + 2; 

a = u(l:n-2); b = u( (n-1) : end) ; 

ap = a; bp = b; 

pos = a(x)max(x,0) ; neg = O(x)min(x,0) ; 
for j = l:length(ii) 
i = abs(ii(j)) ; 
if ii(j) > 
switch(i) 
case 1 

bp(l) = a(l) + pos(b(l)); 
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ap(l) = -b(l) + pos(bp(l)); 
case n-1 

bp(n-2) = a(ii-2) + iieg(b(n-2)) ; 

ap(n-2) = -b(n-2) + neg(bp(n-2) ) ; 
otherwise 

c = a(i-l) - a(i) - pos(b(i)) + neg(b(i-l)); 
ap(l-l) = a(i-l) - pos(b(i-l)) - pos(pos (b(i) ) + c) ; 
bp(i-l) = b(i) + neg(c) ; 

ap(l) = a(i) - neg(b(i)) - neg(neg(b(i-l) ) - c) ; 
bp(l) = b(i-l) - neg(c) ; 
end 

elseif ii(j) < 
switch(i) 
case 1 

bp(l) = -a(l) + pos(b(l)) ; 
ap(l) = b(l) - pos(bp(l)) ; 
case n-1 

bp(n-2) = -a(n-2) + neg(b(n-2)); 
ap(n-2) = b(n-2) - neg(bp(n-2) ) ; 
otherwise 

d = a(i-l) - a(i) + pos(b(i)) - neg(b(i-l)); 
ap(i-l) = a(i-l) + pos(b(i-l)) + pos(pos (b(i) ) - d) ; 
bp(i-l) = b(i) - pos(d) ; 

ap(i) = a(i) + neg(b(i)) + neg(neg(b(i-l) ) + d) ; 
bp(i) = b(i-l) + pos(d) ; 
end 
end 

a = ap; b = bp; 
end 

up = [ap bp] ; 



function int = loopinter(u) 

y.LOQPINTER The number of intersections of a loop with the real axis, 
y, I = LOOPINTER(U) computes the minimum number of intersections of a loop 
7, (encoded in Dynnikov coordinates) with the real axis. (See Moussafir 
7, (2006), Proposition 4.4.) U is either a row-vector, or a matrix of 
7 row-vectors, in which case the function acts vectorially on each row. 

n = size(u,2)/2 + 2; 

a = u(:,l:n-2); b = u( : , (n-1) : end) ; 

cumb = [zeros(size(u, 1) , 1) cumsum(b,2)] ; 

7 The number of intersections before/after the first and last punctures. 
7 See Hall & Yurttas (2009) . 

bO = -max(abs(a) + max(b,0) + cumb( : , 1 : end-1) , [] ,2) ; bn = -bO - sum(b,2); 
int = sum(abs (b) , 2) + sum(abs (a( : , 2 : end) -a( : , 1 : end-1) ) , 2) ... 
+ abs(a(:,l)) + abs (a( : , end) ) + abs(bO) + abs(bn); 



3. proc3_example 

The function proc3_example computes L{u) for a random braid of 4 particles advected 
by the Duffing oscillator, using averaging over 100 realizations (see Section 



function proc3_exainple 
7 Procedure 3 example : 

7 4 particles advected by the Duffing oscillator, average over realizations. 

n = 4; Nreal = 100; 

tmax = 100; dt = 10; npts = 3000; 

X = zeros (npt s, n) ; Y = zeros(npts,n) ; 

LI = zeros (Nreal, tmax/dt) ; tl = dt* (1 : tmax/dt) ; 

rand ( 'twister ' , 2) ; 

for r = 1: Nreal 

f printf (1 ' , 'Realization 7d\n' , r) 



IV Procedure^. 
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FIG. 13: Output of proc3_example (Section A 3). 



% Compute n particle trajectories, from random initial conditions, 
for i = l:n 

[t,xy] = ode45(aduff ing,linspace(0,tmax,npts) ' ,4*rand(l,2)-2) ; 
X(:,i) = xy(:,l); Y(:,i) = xy(:,2); 
end 

[gen,tc] = gencross (t ,X, Y) ; */, Find generators and crossing times. 
'/, Act with the generators on a random initial loop, 
up = randd , 2* (n-2) ) ; up = up/loopinter(up) ; 

for i = 1 : length (gen) , up = [up; loopsigma(gen(i) ,up(end, : ) )] ; end 
y» Find the number of intersections with the real axis. 
L{r} = loopinter (up) ; 

■/, Keep intersections at a list of fixed time intervals dt, for averaging, 
for q = l:(tmax/dt) 

idx = find(tc <= tl(q)); 

Ll(r,q) = L{r}(idx(end)) ; 



logLavg = [0 mean(log(Ll) ,1)] ; tl = [0 tl] ; 
m = polyf it (tl , logLavg, 1) ; 7, Fit a line 
figure(2), plot (tl , logLavg, 'b .-') , hold on 
text (10,1. 6, sprintf ('entropy = .4f ' ,m(l) ) ) 
plot(tl,m(l)*tl+m(2) , 'r-') , hold off 
xlabel('t'), ylabel('<log L>') 

% 

function yt = duffing(t,y) 

delta = 1 ; gamma = 4 ; omega = 2 ; 
yt = zeros (size (y) ) ; 
yt(l,:) = y(2,:); 

yt(2,:) = y(l,:).*(l - y (1 , : ) . *y (1 , : ) ) + gamma*cos (omega*t) - delta*y(2, : ) ; 



end 



end 



